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For active particles the interplay between the self-generated hydrodynamic flow and an external 
shear flow, especially near bounding surfaces, can result in a rich behavior of the particles not easily 
foreseen from the consideration of the active and external driving mechanisms in isolation. For 
instance, under certain conditions, the particles exhibit “rheotaxis,” i.e., they align their direction 
of motion with the plane of shear spanned by the direction of the flow and the normal of the 
bounding surface and move with or against the flow. To date, studies of rheotaxis have focused 
on elongated particles (e.g., spermatozoa), for which rheotaxis can be understood intuitively in 
terms of a “weather vane” mechanism. Here we investigate the possibility that spherical active 
particles, for which the “weather vane” mechanism is excluded due to the symmetry of the shape, 
may nevertheless exhibit rheotaxis. Combining analytical and numerical calculations, we show that, 
for a broad class of spherical active particles, rheotactic behavior may emerge via a mechanism which 
involves “self-trapping” near a hard wall owing to the active propulsion of the particles, combined 
with their rotation, alignment, and “locking” of the direction of motion into the shear plane. In this 
state, the particles move solely up- or downstream at a steady height and orientation. 

PACS numbers: 47.63.Gd, 47.63.mf, 64.75.Xc, 82.70Dd, 47.57.-s 


I. INTRODUCTION 

Flows are ubiquitous in microconfined fluid systems. 
Examples include the flow of blood through capillaries, 
the flow of groundwater through pore spaces in soil, the 
flow of oil through porous rock, and the flow of analyte 
through the channels of a lab-on-a-chip device. Flow 
presents both challenges and opportunities for natural 
and engineered active particles, especially near surfaces. 
For instance, fluid shear can trap bacteria near surfaces, 
inhibiting chemotaxis, but promoting surface attachment 
and biofilm formation pQ. Active particles can harness 
shear flow for sensing and navigation. A motile sper- 
matozoan in shear flow near a surface orients its body 
axis along the flow direction a phenomenon known as 
rheotaxis r] and moves upstream [2] • It is thought that 
in the mammalian oviduct rheotaxis guides sperm to the 
egg 0. Artificial active particles engineered for rheotaxis 
have promising potential applications. For instance, the 
motion of artificial active particles is generally diffusive 
over sufficiently long timescales because thermal noise 
eventually randomizes the orientations of the particles 
4, 5, Rheotactic motion in flow would rectify the di¬ 
rection of motion of the particle. Moreover, by analogy 
with sperm, rheotactic active particles could be used for 
targeted cargo delivery in vivo or in vitro. 


* Corresponding author: uspal@is.mpg.de 

1 In analogy with other “taxis” phenomena, like chemotaxis, 
throughout this paper by rheotaxis we shall refer to the align¬ 
ment of a properly defined director axis of the active particle or 
microorganism along the flow direction. Therefore, depending 
on the relative magnitude of the active motion to that of the 
ambient flow, the net motion of rheotactic particle can yet be 
downstream even if the active component points upstream. 


Experimental studies of the rheotaxis of sperm near 
surfaces have a long history. Bretherton and Rothschild 
investigated live and dead sperm under flow in various 
device geometries 0. They observed that live sperm 
always exhibited motion upstream, i.e., upstream rheo¬ 
taxis. Subsequently, Rothschild observed that motile 
sperm tend to accumulate near surfaces. He attributed 
this to an attraction involving the sperm heads |6]. If 
heads are attracted to surfaces, then reorientation in flow 
can be qualitatively understood as due to the greater 
drag on the tail, which points into the bulk, than on 
the head, which is near the surface. With the head act¬ 
ing as a pivot point, the flow drags the tail downstream, 
like a weather vane 0 0 0. Accumulation of sperm 
and other micro-organisms at surfaces is evidently driven 
by hydrodynamic interactions mm and steric repulsion 
Busmans], although the relative significance of these 
contributions is a subject of debate lanaiEi- 

Recently, the experimental efforts have been extended 
to the study of bacterial rheotaxis. Like sperm, bacte¬ 
ria in quiescent fluid are observed to accumulate near 
surfaces pTOj, nn ua. Hill et al. studied E. coli under 
flow in a microfluidic channel m- They observed that 
the bacteria in the middle of the channel tend to ori¬ 
ent, relative to the flow direction, with a certain acute 
angle. They drift to the sides of the channel, migrating 
across fluid streamlines. Upon reaching the side walls, 
the bacteria orient against the flow and swim upstream. 
As with sperm, these observations could be rationalized 
with a “weather vane” mechanism. More recently, Kaya 
and Koser observed that, within an intermediate band 
of shear rates, E. coli can directly align against the flow 
and swim upstream, without having to migrate to side 
walls [18] . 

Furthermore, we note that for micro-organisms in the 
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bulk (far from bounding surfaces), the interplay of shear, 
swimmer geometry, and swimming strategy can drive ro¬ 
tation of the body axis towards the vorticity direction 
(which is perpendicular to the shear plane), allowing the 
organism to migrate across streamlines. In shear flow, 
bacteria driven by coiled flagella, such as B. subtillis, ex¬ 
perience a lift force in the vorticity direction, owing to 
their chirality. Since this chirality is localized at the flag¬ 
ella, a torque is exerted on the cell body, driving rotation 
towards the vorticity direction [19]. Micro-algae also ro¬ 
tate towards the vorticity direction. The mechanism for 
this orientation is not yet known, but it may be rooted 
in the fact that their two flagella beat at different rates, 
producing unequal propulsive forces and, consequently, a 
torque on the cell body [2D]. This attraction of the ori¬ 
entation vector to the vorticity direction has also been 
labeled rheotaxis [19], but throughout the present study 
we restrict this notion to denote attraction of the orien¬ 
tation vector to the shear plane. 

Theoretical studies have sought to reproduce and iso¬ 
late generic aspects of swimming in flow, including rheo¬ 
taxis. Zottl and Stark considered a model spherical mi¬ 
croswimmer driven by Poiseuille flow in narrow tubes 
and slits [2ll . They found that hydrodynamic inter¬ 
actions with the confining walls can stabilize upstream 
swimming. So-called “pullers” migrated to the center- 
line and swam upstream, while “pushers” were attracted 
to trajectories in which they swim upstream while os¬ 
cillating between the walls. In numerical simulations of 
elongated run-and-tumble particles (posing as model bac¬ 
teria) driven by Poiseuille flow through slit-like channels, 
the particles accumulated near walls due to steric repul¬ 
sion between them and the wall and, due to the subse¬ 
quent alignment of their bodies by the flow, swam up¬ 
stream 122, 23]. Chilukuri et al. modeled swimmers as 
Brownian “dumb-bells” (i.e., two beads connected by a 
spring) |24j . Likewise, when driven by Poiseuille flow 
through a wide slit, a fraction of the particles accumu¬ 
lated near surfaces and swam upstream. This accumula¬ 
tion was enhanced upon including hydrodynamic inter¬ 
actions between the dumb-bells and the walls, indicating 
that, for this system, both steric repulsion and hydrody¬ 
namic interactions promote rheotaxis. Additionally, for 
a Brownian self-propelled particle (i.e., one exposed to 
stochastic forces), ten Hagen et al. found that an un¬ 
bounded shear flow can modify the scaling with time of 
the mean squared displacement of the particle [25]. 

In the last decade, significant progress has been made 
in developing synthetic self-propelled particles. These 
particles are envisioned as key components in future lab- 
on-a-chip, drug delivery, and pollution remediation sys¬ 
tems [261 El]- In these and other applications, active 
particles will undoubtedly encounter shear flows near sur¬ 
faces and will need to autonomously sense and respond to 
flow. Many generic aspects of artificial microswimmers 
in flow are likely captured by the theoretical and numer¬ 
ical studies discussed above. However, achieving rheo¬ 
taxis in artificial systems, for which active stresses can 


be generated through, e.g., local changes in the chem¬ 
istry of the environment, requires an in-depth analysis 
of the interplay between the external flow, the presence 
of bounding surfaces, and the mechanism of propulsion 
(i.e., activity). For example, for a large class of active 
particles, a reaction involving a “fuel” molecule in the 
surrounding solution is catalyzed by a region of the par¬ 
ticle surface, leading to propulsion either via the forma¬ 
tion and expulsion of bubbles [21], or via “self-phoresis,” 
i.e., the generation of a tangential surface flow powered 
by the gradient of the product and/or reactant molecules 
[221 EDI- For the latter case, we have recently shown that 
the complex behavior emerging when the particle motion 
occurs near confining boundaries can be understood only 
if the coupling loop between the distributions of chemical 
species and hydrodynamics is fully accounted for |3T . 

There have been comparatively few studies of artifi¬ 
cial microswimmers in flow, and the issue of rheotactic 
behavior was rarely raised. Sanchez et al. have shown 
that bubble-powered, rod-shaped catalytic “microjets” 
can move upstream in a channel when guided by an ex¬ 
ternal magnetic field [21]. Frankel and Khair studied 
theoretically a model of self-phoretic Janus particles and 
have shown that drift across streamlines in unbounded 
shear flow may occur at finite Peclet number [E2J- Tao 
and Kapral have used numerical simulations in order to 
investigate the motion of a self-phoretic dimeric nanomo¬ 
tor moving upstream against a Poiseuille flow in a chan¬ 
nel and to examine the effect of confined flow on the 
propulsion velocity of the motor [33] . However, they 
used a potential which confined the motor to the channel 
centerline, and therefore did not probe the possibility of 
rheotaxis. Recently, Palacci et al. employed microflu¬ 
idic chips and self-propelled dimers made out of polymer 
spheres and hematite cubes partially embedded in the 
polymer matrix in order to provide the first experimen¬ 
tal demonstration of rheotaxis, via the “weather vane” 
mechanism, for an artificial swimmer [34] . The heavy 
dimers sediment to the bottom surface and, once being 
exposed to blue light and to hydrogen peroxide “fuel”, 
the hematite becomes catalytically active and promotes 
the decomposition of the hydrogen peroxide. This leads 
to an attraction of the hematite end to the surface and 
thus, while the whole dimer tends to move due to the 
concentration gradients, the hematite end serves as an 
“anchor” for the dimer to orient with the ambient, exter¬ 
nally controlled flow within the microfluidic chip. 

In the present study we address the issue of whether 
artificial spherical active particles can exhibit rheotaxis 
when moving in the vicinity of a planar surface. Without 
a distinction between major and minor axes, the particle 
cannot align via the intuitive “weather vane” effect. Thus 
it is not a priori clear if rheotaxis is possible for such 
systems; for rheotatic behavior to emerge nonetheless a 
different physical mechanism must be at work. Consid¬ 
ering the overwhelming importance of spherical particles 
for experiments and applications of artificial swimmers 
(a spherical particle is naturally simpler and easier to 
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FIG. 1. (a) A spherical active particle is suspended in a linear shear flow above a planar wall which is translationally invariant 

in x and y directions. Here, it is shown as a Janus particle which moves away from its catalytically active cap (blue). The 
orientation of the particle is indicated by the unit vector p (green), which is described by the angles 9 and <j>. Black arrows 
and lines correspond to projections onto the xy and x'y' plane, (b) Illustration of the “primed” coordinate frame used to solve 
subproblem II (see main text). In both (a) and (b) the planar wall is located at z = 0 and z' = 0, respectively. 


fabricate than an elongated body), this question out of 
basic research turns out to be one of significant practical 
interest, too. Furthermore, the answer to this question 
may also shed light on the motility in confined flows of 
micro-organisms with spherical or nearly spherical shape. 

We focus our theoretical and numerical study on the 
case of spherical particles with an axisymmetric propul¬ 
sion mechanism which can be described via an effective 
surface slip velocity. This framework captures a number 
of important classes of microswimmers, including self- 
phoretic catalytically active Janus particles, as well as the 
classical “squirmer” model, introduced by Lighthill [35] 
and subsequently refined by Blake >55] . which is known 
to account very well for the essential features of the self¬ 
propulsion of ciliated micro-organisms. We show that if 
the swimming activity of the particle can lead to a “self¬ 
trapping” near a bounding surface rheotactic behavior 
may emerge. We establish the conditions under which 
this “self-trapping” is accompanied by a mechanism for 
resisting rotation by the flow and maintaining a steady 
orientation. In contrast with the “weather vane” mech¬ 
anism, rotation of the orientation vector (i.e., the sym¬ 
metry axis of the particle) into the shear plane is now 
driven not by the external flow, but rather by the near¬ 
surface swimming activity of the particle. Building on 
these results, as well as on our recent study of the self¬ 
propulsion of a catalytic Janus particle near a surface in 
a quiescent fluid GU, we show how the surface chemistry 
of a catalytically active Janus particle can be designed to 
achieve rheotaxis. The wide applicability of our results is 
further emphasized by demonstrating that rheotactic be¬ 
havior also occurs for squirmers which in quiescent fluids 
are attracted by surfaces mm- 


II. THEORY 

A. Problem formulation 

We consider a spherical, axisymmetric, active particle 
suspended in a linear shear flow in the vicinity of a pla¬ 
nar wall (Fig. 0a)). We adopt a coordinate frame (the 
“lab frame”) in which the wall is stationary. The wall is 
located at 2 = 0 and has a normal vector in the z direc¬ 
tion. The particle has radius R and its center is located 
at r p = (x p ,y p , z p ) T ; h = z p denotes the height of the 
particle center above the wall. We describe the particle 
orientation in terms of the “director” p, |p| = 1, which 
is the direction of active motion if the particle would be 
suspended in an unconfined, quiescent fluid. Due to ax¬ 
ial symmetry, p is parallel to the axis of symmetry of 
the particle. (For instance, for a catalytic Janus particle 
which moves “away” from its catalytic cap, the direction 
of p is given by the vector from the pole of the catalytic 
cap to the pole of the inert region.) The height and 
orientation completely specify the instantaneous particle 
configuration. We shall find it convenient to alternatively 
describe the particle orientation in terms of spherical co¬ 
ordinates 9 and (/>, where 6 is the angle between p and 
z, and (j> is the angle between x and the projection of p 
onto the xy plane (see Fig. 0a)). We seek to calculate 
the translational velocity U and angular velocity f 1 of 
the particle as a function of the configuration (h, p). 

The active particles we consider here are such that 
the self-propulsion mechanism can be accounted for via 
an effective surface slip velocity v, s . The slip velocity 
can encode, for example, the surface flow generated by a 
“squirming” mechanical deformation of the particle sur- 
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face (e.g., the time averaged motion of the cilia of a micro¬ 
organism), or the surface flow generated by the interac¬ 
tions between reactant and product molecules with the 
particle (as in the case of a catalytically active Janus 
colloid). Here we take the slip velocity as a known func¬ 
tion of the position along the particle surface and of the 
configuration (ft, p). We postpone for Section II.E the 
discussion of how to obtain this slip velocity for a given 
swimmer model. 

The velocity in the suspending fluid is u(r) = u ea;t (r)-|- 
u d (r), where u ext (r) = jzy is the external background 
flow, Ud is the disturbance flow created by the parti¬ 
cle, and 7 denotes the spatially and temporally constant 
shear rate. The fluid velocity is governed by the Stokes 
equation —VP + rjS7 2 \i = 0, where P(r) is the pressure 
field and p is the viscosity of the solution, as well as the 
incompressibility condition V ■ u = 0. On the planar 
wall, the velocity satisfies the no-slip boundary condition 
u = 0. On the surface of the particle, accounting for 
the effective slip implies that the fluid velocity satisfies 
u = U + S7x (r — r. p ) + v s . Finally, we specify that the 
active particle is force and torque free, thus closing the 
system of equations for the six unknowns U and Cl. 


of h/R in Ref. [55] , 

We now turn to subproblem II. Due to the symme¬ 
try of the system and the absence of thermal fluctua¬ 
tions, the motion of the particle is confined to the plane 
spanned by the orientation vector p and the wall normal 
z. Therefore it is convenient to introduce a “primed” 
reference frame to solve this subproblem (Fig. [ljb)) . In 
the primed frame, the wall is stationary and located at 
z' = 0. The vector z' is normal to the wall, the vec¬ 
tor y' lies in the plane containing the particle symmetry 
axis and z' (which means that the projection of p onto 
the x'y' plane points into the direction of y'), and x! 
is orthogonal to y' and z!. In general, this subproblem 
must be solved numerically. We discuss the numerical 
solution in Section III. Here we assume the solution to 
be known and we use the most general form compatible 
with the symmetry constraints discussed above, i.e., a 
translational velocity = U a , y 'y' + U a , z 'z' and an an¬ 
gular velocity Cl' a = Cl a yx', where U a>y >, U a ^, and Cl a ^ 
are functions of h/R and 9. Transforming back to the 
original reference frame, we obtain 

U a = TJa.y’ COS {(j))x + U a y sm((/))y - 1 - U a ,z'Z , 

Cla = Cl a ^ sm((/))x - Cl a ^ cos {cj))y . (2) 


B. Solution by linear superposition 

We exploit the linearity of the Stokes equation in order 
to split the complete problem into two subproblems. In 
subproblem I we consider the motion of an inert (non¬ 
active) sphere with configuration (ft, p) driven by an ex¬ 
ternal shear flow. We note that although, for consistency 
of notations, p is included as a configuration variable, in 
the case considered in this subproblem the velocity of the 
particle actually does not depend on p (but, obviously, 
p will evolve in time due to rotation of the particle.) In 
subproblem II, a spherical active particle with instan¬ 
taneous configuration (ft., p) moves through a quiescent 
fluid. The particle velocity for the complete problem is 
then obtained by linear superposition: U = U f + U a and 
f2 = Clf + Cl a , where the velocities U / and Cl f are due to 
/low (subproblem I), and the velocities U a and Cl a are 
due to the activity-induced motion (subproblem II). 

Subproblem I was solved analytically by Goldman et 
al. using bispherical coordinates [32]. The translational 
and rotational velocities of the particle are given by 

U / = 7 hg(h/R)y, 

Cl f = -\if{h/R)x. (1) 

Note that the external flow does not cause motion of 
the particle in the z direction. The functions f{h/R ) 
and g(h/R) encode the hydrodynamic interaction with 
the no-slip wall, which retards both the translation and 
the rotation of the particle. These functions increase 
monotonically with h/R and have the properties 0 < 
f(h/R) < 1, 0 < g(h/R) < 1, f{h/R —> 00) —> 1, and 
g{h/R —» 00) —»• 1. They are tabulated for selected values 


C. Particle motion and fixed points for the 
dynamics 


In the absence of thermal fluctuations (which we as¬ 
sume to be negligible), the time evolution of the particle 
configuration (ft, p) is determined, in the overdamped 
limit, by the translational and rotational velocities de¬ 
rived above. By noting that the motion of the director 
obeys p = Cl x p and that (see Fig. [TJa)) 


sin(</>) 


Py 

V 1 pI ’ 


cos(^) 


Px 

X / 1 -Pz ’ 


we arrive at the following dynamical equations (using 
fl = flf + fi a , U = TJf + U a , p = ft x p): 


Px — QyPz ^zPy — 

Py = ^zPx ^xPz 


^ a,x' (jPz, hjR)p x Pz 

V l ~P 2 z 


1 . ^a,x'(Pzi h/R)p y p z 

= if1Pzf{h/R) - - 


V 1 ~P 2 z 


Pz — ^xPy ^yPx 


and 


= -^iPyfih/R) + tt a ,x'(Pz, h/R)y/l -p 2 z 


h =u z = U a ^{p Zl h/R ), 


(3) 

(4) 

(5) 

( 6 ) 


where the dot over a quantity denotes its time derivative. 
Since p 2 +p 2 +pl = 1, the dynamics of the director must 
obey the constraint p x p x PyP y + p z Pz = 0; this is ob¬ 
viously satisfied by Eqs. ([3])-|5]). Since, accordingly, the 
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three components of p are not independent, the dynam¬ 
ical system is de facto three-dimensional. 

Defining the generalized velocity vector V = 
[Px,Py,Pz, h] T , we search for fixed points (p *,h*) of the 
dynamics, which are determined by V(p*, ft*) = 0. Such 
configurations (p *,h*) correspond to the particle trans¬ 
lating along the wall with a fixed height and orientation, 
i.e., steady states with properties which are compatible 
with (although not all of them are necessary for) rheo- 
taxis. In order to obtain (p* ,h*), we start by considering 
the three possibilities for p x = 0 to be satisfied. For com¬ 
pleteness, we discuss here all three cases and, for each 
of them, derive the additional conditions following from 
p y = 0, p z = 0, and ft. = 0. In anticipation, we note that 
for the types of swimmer considered here, it turns out 
that only the third case, which corresponds to p* = 0 , is 
compatible with rheotaxis. Briefly, in the first two cases, 
the combined effects of shear and activity produce a fixed 
point with p* = 0, i.e., 9 = 90°. Most types of swimmer 
do not satisfy the restrictive conditions required to pro¬ 
duce such a fixed point. Consequently, the rest of the 
study in this paper will focus on the third case. 

(i) The case f l a ,x'{p z = p 2 ,ft*) = 0. This is a particle 
configuration for which in the absence of external flow the 
particle only translates. In order to also satisfy p y = 0, 
p z = 0 , and ft, = 0 , it is then required that p* = 0 , p* = 0 , 
and U a z ' (p* = 0, ft*) = 0. The first two of the latter re¬ 
lations imply that the director is along the x-axis, which 
coincides with the vorticity axis of the external shear 
fiov 0 thus the director is not rotated by the external 
flow, while the entire particle rotates around the director 
due to the external shear flow, which also advects the par¬ 
ticle (a “logrolling” state). Furthermore, withp* = 0 the 
conditions il' a (p* = 0,ft*) = 0 and U a , z '{pt = 0, ft*) = 0 
imply that there must be a fixed point for subproblem II 
which occurs at p* = 0, i.e., at 9* = 90°, for which the 
director is parallel to the wall (see Fig. 0 a)). We note 
that this condition on subproblem II is very restrictive; it 
turns out that it cannot be satisfied by any of the active 
particles considered in Sec. III. 

(ii) The case p* = 0. In this case the director lies 
within a plane parallel to the wall (see Fig. 0a)). The 
equation p y = 0 is then automatically satisfied (Eq. <@»- 
The condition ft = 0 implies that the curve p z (ft) defined 
by U atZ '(p z , ft) = 0 passes through 9 = 90° (i.e., p z = 0) 
at a certain height ft*. The condition p 2 = 0 then implies 

p* = ~ ■ The corresponding physical picture 

is as follows. As the director is within a plane parallel 
to the wall, the contributions to p from both shear and 
activity are entirely along the z direction: for p z = 0 
Eqs. (0 - (0 imply p = (0, 0,p 2 ). The magnitude of 
the contribution from the shear flow (~ 7 in Eq. ©) 


2 The vorticity u of the flow u is defined as u = V X u; for the 
shear flow u ex t = fzy one has uj = - A rx so that the vorticity 
direction is thus everywhere aligned with the ai-axis. 


depends on the angle <j>. For instance, if the director is 
oriented along the x-axis (p y = 0 , so that <j> = 0, tt), the 
contribution of shear to p 2 is zero (see Eq. ([5])). At p y = 
p* (see above, i.e., for a certain angle cj) 7 ^ 0 , tt), due to the 
definition of p* the contributions from shear and activity 
to Pz sum to zero. For this fixed point to occur, due to 
|p| = 1 it is required that |p*| = 2n a.*'( p *-°’ h 1 


< 1 . 


•yf(h*/R) 

This requirement, together with that of the vanishing of 
U a , z ' at a certain height ft* for 9 = 90°, places weaker 
demands on subproblem II than the requirements in case 
(i). Nevertheless, as in case (i), it turns out that none 
of the types of active particles which will be considered 
in Sec. Ill satisfies the fixed point requirement on I/ a z /, 
and thus none of them can fit case (ii) either. 

(iii) The case p* = 0. In this case the director p lies 
within the shear plane (i.e., the yz plane in Fig. 0 a)). 
If a fixed point exists, in that state the particle trans¬ 
lates either only upstream or only downstream because 
once having reached the fixed point of the dynamics the 
orientation cannot vary in time. If the fixed point is sta¬ 
ble (this issue will be addressed analytically in the next 
subsection and numerically in Sec. Ill), this state corre¬ 
sponds to rheotaxis. Since in this case p y = ± \J\ — p 2 
and thus p y p y = —p 2 p 2 , Eqs. 0 and (0 are identi¬ 
cal (as expected, because only two of the Eqs. 
are independent). Therefore, the steady height ft* and 
orientation p* are obtained as the solutions of 


sgn (p* y ) n a y(p* z ,h*/R) = -if(h*/R) 
U a<z ,(p* z ,h*/R) = 0. 


(7) 


The first equation simply expresses that at the steady 
state the sum of the contributions to the angular veloc¬ 
ities from shear and activity must add up to zero. The 
term sgn(p*) appears because shear tends to rotate the 
director away from the wall when p y < 0 (i.e., <f> = 270°) 
and towards the wall when p y > 0 (i.e., <f> = 90°) (see 
Fig. 0a) and recall that in the present case p lies in 
the yz plane), whereas, owing to the symmetry of sub¬ 
problem II, there is no dependence on <p for the direction 
of rotation by the swimming activity. In contrast to the 
cases (i) and (ii), it turns out that the conditions in Eqs. 
0 can be satisfied by all three classes of spherical ac¬ 
tive particles considered here (see Sec. III). Therefore we 
proceed further with the analysis of the case p* = 0 . 


D. Linear stability analysis of the fixed point p* = 0 


In order to establish the conditions under which the 
fixed point p* = 0 (case (iii) above) is stable, so that the 
corresponding steady state exhibits rheotaxis, we proceed 
by carrying out a linear stability analysis. Since for p x = 
0 the director lies within the shear plane and, due to 
the symmetry of the problem, cannot rotate out of this 
plane, in the three-dimensional phase space defined by 
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FIG. 2. (a) Schematic description of the particle dynamics if the director p = ( Px,Py,Pz ) (green) lies in the indicated shear 

plane. The angle £ is formed by the director and the flow direction. The wall is located at 2 = 0 and h is the distance of the 
particle center from the wall. For a catalytic Janus particle, the angle ip characterizes the size of the catalytic cap (blue). Note 
that here ip is larger than in Fig. [I] (b) Illustration of a particle in a rheotactic steady state which occurs for p* < 0 and 
pi < 0, i.e., for £ within the range n < £ < 37 t/ 2 . The external flow contributes a clockwise component — ^ 7 f(h/R)x (cyan) 
to the particle angular velocity (Eq. Q). For the angle £ to remain stationary, the swimming activity of the particle must 
contribute a counterclockwise component (magenta) with equal magnitude to the angular velocity. Since the particle, assumed 
to be axisymmetric, does not rotate in free space, i.e., without the wall and in the absence of the external shear flow, this 
contribution must be due to the interaction with the wall. If the swimming velocity of the particle is large enough, the particle 
swims upstream, (c) Upstream rheotaxis is stable against perturbations of the director out of the shear plane, where p x = 0. 
Shear flow tends to rotate the director around the x direction, keeping p x constant and thus does not contribute to p x (see the 
last paragraph in Subsec. II.C). The swimming activity of the particle tends to rotate the orientation vector along “lines of 
longitude” on the unit sphere |p| = 1. For a small perturbation A p x (blue line segment) away from the rheotactic state shown 
in (b), activity rotates the director (green vector) towards the pole p z = —1, decreasing p x and p y due to |p| = 1. The effect 
of activity is indicated by the magenta arrow. 


p x , Pz, and h a trajectory with initial condition in the 
plane p x = 0 will be confined to that plane for all times 
t. Upon linearization of the equations of motion near the 
fixed point, a small perturbation A p x decouples from the 
other variables: 


dAp x _ n a ^(p* z ,h*/R)p* z 
dt -pf 


A p x . 


After using Eq. |7]) for substituting 

sgn (pl)n a ^(p* z ,h*/R) = 57 f(h*/R), we obtain 


dAp x jf(h*/R)p* z 


dt 


2 pi 


Ap x . 


Therefore Ap x decays or grows exponentially, 

-7 f{h*/R)pl 


A p x = Ap x (0) exp 


2 Pi 


( 8 ) 


(9) 


( 10 ) 


depending on the sign of the quantity multiplying t in 
the argument of the exponential. 

Since f(h*/R ) > 0 and 7 > 0, the stability of the fixed 
point is determined by the sign of the ratio p* z /p y ■ For 
the particle to swim oriented upstream (i.e., to exhibit 
“positive” rheotaxis), it must align against the flow, so 
that p^ < 0 ; stability then requires that p* z < 0 , i.e., 


the director points towards the wall. Additionally, if the 
swimming velocity is large enough to overcome the shear 
flow, i.e., for y p < 0, where y p = —U ay i + jh*g(h*/R) 
(see Uf. y in Eq. ([T]) and U a , y in Eq. ([ 2 j with cp = 270°), 
net upstream swimming emerges. The other possibility, 
that an attractor exists for which p* and p* z are both pos¬ 
itive, is in general unlikely, because the particle would be 
pointing away from the wall (see Fig. |TJa) ). Further¬ 
more, it would correspond to the less interesting case of 
negative rheotaxis, in which the particle motion is in the 
flow direction. Therefore, it will not be discussed further. 


The analysis of the linear stability shows that we need 
to examine particle dynamics only in the symmetry plane 
p x = 0 in order to determine whether rheotaxis occurs. 
In the following, we shall examine the two-dimensional 
dynamics in this plane in detail. Additionally, we shall 
perform fully three-dimensional numerical calculations of 
trajectories with initial conditions p x (t = 0 ) yf 0 in order 
to explore the limits of validity of the linear analysis, as 
well as to understand the details of the evolution towards 
alignment in the shear plane. We introduce a new angle £, 
shown in Fig. |2ja) , which will turn out to be convenient 
for describing the director orientation if p x = 0. The 
requirement that the attractor occurs with p* < 0 and 
p* z < 0 translates into n < < 3n/2. 

Before proceeding with the analysis, we discuss the 
physical picture behind the stability criteria in Eq. (10 1 . 
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While we have provided mathematical conditions for 
rheotaxis, the corresponding physical interpretation can 
provide further insight. Equation <§ shows that particle 
activity has an effect on p x , but shear does not. The con¬ 
tribution of shear to the angular velocity of the particle is 
rigid rotation around the x axis. In the absence of activ¬ 
ity, shear rotates the tip of the director in a coaxial circle 
around the x axis, keeping p x constant. Now we consider 
activity. Due to the mirror symmetry with respect to the 
plane containing p, the contribution of activity can only 
drive rotation of the director with arbitrary orientation p 
towards or away from the planar wall. On the spherical 
surface defined by |p| = 1 , this corresponds to rotation 
of the director along a line of longitude, where p z = 1 
and p z = —1 define the poles [Fig. |2jc)]. The question, 
then, is whether, for a small perturbation A p x , activity 
will tend to drive the director towards the equator p z = 0 
of |p| = 1 or towards the nearest pole. From Fig. |2j[c) 
we note that the effect of the small perturbation A p x is 
to shift the tip of the director to a neighboring line of 
longitude. Therefore, if activity tends to rotate the di¬ 
rector towards the equator, A p x will increase (see Fig. 
He)) because the distance between neighboring lines of 
longitude increases as they approach the equator from 
the pole. The direction of this activity induced rotation 
(i.e., towards the pole or towards the equator) is deter¬ 
mined by the value of f *. For instance, consider a positive 
rheotactic state, for which n < £* < 3n/2. In order to 
balance rotation by shear, activity must tend to rotate 
the director towards the pole p z = —1, as shown in Fig. 
[2jb) and (c). Hence, this dynamical fixed point is stable 
against small perturbations in p x . 


E. Calculation of the slip velocity 

Finally, we briefly outline the calculation (if necessary) 
of the slip velocity v s = v s ^ p e' e ^ (where e' g denotes the 
unit vector in the primed coordinate frame corresponding 
to the polar angle 9 p measured from the director, i.e., the 
symmetry axis of the particle) along the surface of the 
spherical particle which, for the specific swimmer models 
considered here, is an input to the problem outlined in 
Subsec. II.A. 

For a “squirmer”, v s is specified by fiat and does not 
depend on the configuration (p, h) of the particle. Fol¬ 
lowing Li and Ardekani [35], we consider a squirmer for 
which the slip velocity is given by 

v s ,B p = Bi sin (6 p ) + B 2 sm(d p ) cos (9 P ), (11) 

The model is characterized by the squirming mode am¬ 
plitudes B i and B 2 . In an unconfined quiescent fluid 
(i.e., free space) such a squirmer moves with velocity 
U f .,. = (2/3) A?!. 

A self-phoretic Janus particle generates solute 
molecules over the surface of a catalytic cap. We param¬ 
eterize the extent of the cap by \o = — cosfip), where 


the angle if is defined in Fig. Ha). It is assumed that 
the interaction between the active particle and the (non- 
uniformly distributed) solute molecules has a range much 
smaller than R and h. Thus it can be accounted for as 
driving a tangential flow in a thin boundary layer of thick¬ 
ness 5 « R,h surrounding the particle. This can be 
modeled as an effective slip velocity v s = — 6 (r)V||c(r), 
where c(r) is the solute number density field (2j2 [3D] • The 
operator V|| = (1 — nn) • V is the projection of the gra¬ 
dient onto the active particle surface, with n denoting 
the surface normal oriented into the fluid. The quantity 
b{ r) is the so-called “surface mobility,” determined by 
the details of the interaction between the active parti¬ 
cle and the solute molecules. For instance, the sign of b 
expresses whether the active particle is attracted to or re¬ 
pelled from the solute molecules. Throughout the present 
study, we assume repulsion. Since the solute number den¬ 
sity c(r) is affected by the presence of the wall, in this 
model the slip velocity v s depends on the configuration 

(p,/i). 

If the effect of advection on c(r) can be neglected, one 
can determine c(r) (and with this v s ) prior to the con¬ 
sideration of the hydrodynamic problem posed in Subsec. 
II.A (see, e.g., Refs. [GD/SO])- (Within this approach it is 
disregarded how the flow field u couples back to c(r).) In 
this case, the solute number density is quasi-static, and 
obeys the Laplace equation DV 2 c = 0, subject to the 
boundary conditions for the normal derivatives h- Vc = 0 
on the wall, h ■ Vc = 0 over the inert surface of the col¬ 
loid, and — D{n-Vc ) = n on the catalytic cap. The quan¬ 
tity D is the diffusion coefficient of the solute molecules, 
and k is the flux of solute per unit area from the cat- 
alytically active cap. Concerning the first and second 
boundary conditions, the wall and the inert surface of 
the particle are taken to be completely impenetrable to 
solute molecules. In the third boundary condition, it is 
assumed that the rate of the reaction is independent of 
the local concentration of solute molecules. 

It is valid to neglect advective effects if the Peclet 
numbers Pe = Uf. s .R/D and Pe y = fiR 2 /D are small, 
where Uf, s _ is the particle velocity in free space. (The 
Peclet number Pe is the ratio between the time scale 
for solute to diffuse a distance R and the timescale 
for the particle to move a distance R. Similarly, the 
Peclet number Pe-y characterizes the advection of the so¬ 
lute molecules by the external flow, for which the ve¬ 
locity scale is yi?, relative to the diffusion of the so¬ 
lute molecules [32]•) Additionally, we assume that the 
Reynolds number Re = pUf. s .R/p (where p is the mass 
density of the suspending fluid and p is its dynamic vis¬ 
cosity) is small, permitting the use of the Stokes equa¬ 
tions to account for the hydrodynamics. For example, a 
self-propelled Janus colloid which catalyzes the decompo¬ 
sition of hydrogen peroxide into water and oxygen typ¬ 
ically has a diameter of ~ 5 pm and moves with a ve¬ 
locity ~ Uf.s. = 5 pm/s through the aqueous solution 
(p ~ 10 3 kg/m 3 and p ~ 10 -3 Pa x s). In this case, 
Re ~ 10 -5 . Furthermore, estimating the diffusion coefR- 


dent of oxygen in water as D ~ 4 x 10 -9 m 2 /s, we obtain 
Pe ps 10 -3 [41]. A representative shear rate in microflu¬ 
idic devices is 7 = 1 s -1 (see, for instance, Refs. [5] and 
[20]). which renders Pe-y ss 10 -3 . Thus for typical active 
particles the above assumptions of small Pe, Pe -y, and 
Re numbers are well justified. 

Before proceeding to an analysis of the numerical re¬ 
sults, we remark on the physical realizations which would 
be captured by our model. In many experiments, the cat¬ 
alytic reactions involve more than one product as well as 
possibly a number of reactants. If the system remains in 
a reaction-rate limited regime (i.e., the reactants are in 
abundance and transported sufficiently fast so that there 
is no noticeable depletion near the catalyst), accounting 
for more than one product means to modify the expres¬ 
sion for the slip velocity in a linear manner: the gradient 
of each reaction product “i” multiplied by its correspond¬ 
ing surface mobility by is simply added in order to obtain 
the phoretic slip around the particle. In this case, the 
results presented here can be easily extended by replac¬ 
ing the effective number density c by the sum over the 
densities of all products with an effective b such that 
6 V||C —► foiV 11 ci + & 2 V||C 2 + .... On the other hand, if 
the catalytic reaction is diffusion limited and involves at 
least two reactants, the source boundary condition at the 
catalytic-active region involves products of the densities 
of the reactants. In this case the equations describing 
the system become nonlinear, and the present model is 
neither applicable, nor does it lend itself to an obvious 
extension. 


III. NUMERICAL RESULTS AND DISCUSSION 


When considering a Janus particle, we use the bound¬ 
ary element method (BEM) in order to solve the diffusion 
equation for the solute number density field c(r). Addi¬ 
tionally, for both propulsion mechanisms (the squirmer 
and the catalytically active particle), we use the BEM 
in order to solve the hydrodynamic subproblems I and 
II. A detailed introduction to the BEM is provided by 
Pozrikidis [42] . We adopted his freely available BEMLIB 
code and modified it for the present study. For subprob¬ 
lem I, we calculate f(h/R) and g(h/R) for an (inactive) 
sphere in a shear flow at various heights h/R. We find 
that the numerically obtained values are in good agree¬ 
ment with the analytically obtained values given by Gold¬ 
man et al [231 ■ In a previous study, we provided a de¬ 
tailed description of how to apply the BEM in order to 
solve the diffusion equation (for a catalytic Janus par¬ 
ticle) and subproblem II [31]. Within this approach we 
have obtained U a , z ' and f l a , x i on a grid of h and 9. 

In order to determine full particle trajectories for an 
initial condition (ho, po), we interpolate f(h/R ), g(h/R), 
U a , z ', and f 1 a ,x' and integrate Eqs. <§-©> numerically. 
In the following, we restrict our numerical calculations 
to the region 1.02 < h/R < 7.1. Above h/R = 7.1, we 
consider the particle to have “escaped” the surface. Be- 



“hoverer” 

large catalytic cap (blue) 

achieves motionless “hovering” states near 
surfaces, with the cap oriented vertically away 
from the surface 

“slider” 

half coverage by catalyst 

different strengths of interaction of the solute with 

the catalytic (blue) and inert (red) “faces” 

“slides” along surfaces at a steady height and 
orientation 



“squirmer” 

slip velocity specified by fiat as a combination of 
squirming modes 

can “slide” near surfaces for certain combinations 
of squirming mode amplitudes 

model for the swimming of ciliated micro-organisms 


FIG. 3. Illustration of the three types of spherical active par¬ 
ticles studied here numerically. The three types differ in the 
kind of their surface activity. A catalytically active particle 
generates solute molecules (green discs) over the surface of a 
catalytic cap (blue). The interaction between the solute and 
the particle surface, here taken to be repulsive, drives a sur¬ 
face flow (slip velocity) as discussed in the main text, leading 
to a motion of the particle away from its cap. A “hoverer” 
(Sec. |III A I exhibits a very high coverage by catalyst. As 
discussed in Ref. eh, if it is near a planar wall it can achieve 
a motionless “hovering” state at a steady height above the 
wall and with the cap oriented away from the wall (as shown 
in Fig. [ 4 ] left panel). For the hoverer considered here, we 
take the solute to interact identically with the catalytic and 
inert “faces” of the particle. In contrast, the “slider” consid¬ 
ered in Sec. EH is only half covered by catalyst, but the 
solute interacts more strongly with its catalytic face (blue) 
than with its inert face (red). This particle tends to swim 
along surfaces at a steady height and orientation, or “slides.” 
Finally, for the “squirmer” (Sec. |III C I, the slip velocity is not 
due to a self-generated distribution of solute, but rather it 
is specified by fiat as a combination of “squirming modes.” 
For certain combinations of squirming mode amplitudes, the 
squirmer can achieve sliding state, too. 


low h/R = 1.02, various of our physical approximations 
(e.g., that the effects of the solute distribution around 
the particle can be accounted for by a phoretic slip cal¬ 
culated within a thin layer, including that part of the 
surface of the particle which is in close proximity of the 
wall) are expected to break down. For a catalytic Janus 
particle, Uq = \b\n/D provides a characteristic velocity, 
and T 0 = R/Uq a characteristic timescale. For instance, 
it was shown analytically that a catalytic Janus particle 
with uniform surface mobility and half coverage by cat¬ 
alyst has a velocity Uf. s ./Uo = 1/4 in free space [44]. 
When we take the catalytic cap and the inert region to 
have unequal surface mobilities, we use Uo = \b cap \n/D. 
Similarly, for squirmers the amplitude Bi of the first 
squirming mode provides a characteristic velocity scale, 
and the time T s = R/By a characteristic timescale. The 
ratio between the squirming mode amplitudes defines a 
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FIG. 4. Schematic illustration of a dynamical process through 
which a “hoverer” particle can achieve a rheotactic state. Left 
panel: The particle is initially in a “hovering” state without 
any external flow. As discussed in Ref. ED. the catalytic 
cap of the particle (blue) generates solute molecules (green 
discs), which interact with the particle surface via a repulsive 
potential. A “cushion” of repulsive solute in the space be¬ 
tween the particle and the hard planar wall exactly balances 
the tendency of the particle to swim away from its cap. Cen¬ 
ter panel: An external shear flow tends to rotate the particle 
clockwise (cyan arrow). The near-wall chemical activity of 
the particle tends to rotate the particle back to the hovering 
orientation (magneta arrow). However, this contribution to 
particle rotation is small for a small displacement from the 
hovering orientation. The net effect is that the particle ro¬ 
tates away from the hovering orientation. Right panel: For 
sufficiently large displacements of the angle £ from the hov¬ 
ering state value £ = 37r/2, the two contributions to particle 
rotation balance, and a steady angle is achieved. Simultane¬ 
ously, the particle still achieves a steady height above the wall 
due to the solute “cushion” effect. As a result, the particle 
translates upstream (black arrow), as its catalytic cap (blue) 
is oriented slightly downstream. 


parameter B 2 /B 1 . In each of the cases which we shall 
study in this section, the corresponding velocity and time 
scales defined above and the particle radius R will be em¬ 
ployed to render the dynamical equations, Eqs. 
dimensionless. 

First, we shall apply the theoretical results derived in 
the previous section to show that the surface chemistry of 
a catalytically active Janus particle can be tailored such 
that it leads to the occurrence of positive (upstream) 
rheotaxis. We shall provide two rather distinct examples 
of such a design of the surface properties, each exploit¬ 
ing a particular pathway to produce the stabilizing wall- 
induced rotation component discussed in Subsec. |ITd| 
These two examples are depicted schematically in Fig. 
[3] In order to illustrate the generality of our theoretical 
results, we shall show that rheotaxis can occur for certain 
spherical squirmers. 

Our approach to design is guided by the idea outlined 
in Fig. [2jb) . For positive rheotaxis, the particle director 
p must point upstream and towards the wall. The exter¬ 
nal flow contributes a clockwise component, shown by the 
cyan arrow, to the angular velocity O of the particle. For 
the particle to maintain a steady orientation, the near¬ 
wall swimming activity of the particle must contribute a 
counterclockwise component fl a = — Clf, shown by the 
magenta arrow, to f2 = rij + fi a . Since the axisymmetric 
particle does not rotate in free space, the counterclock¬ 


wise component Cl a must be due to the effect of the wall 
on the fluid velocity field and the solute number den¬ 
sity field. Additionally, the particle must be attracted 
to a steady height through its near-wall swimming ac¬ 
tivity. In other words, the tendency of the particle to 
swim away from its cap must be counteracted by the ef¬ 
fect of the wall on the fluid velocity and solute density 
field. In the following two subsections, we shall introduce 
two particle surface chemistries which fulfill these crite¬ 
ria. By tailoring the surface chemistry, we can turn on or 
off, and rationally control, various physical mechanisms 
which contribute to the particle motion. 

We note that the following three subsections (corre¬ 
sponding to the three particle designs shown in Fig. [ 3 ]) 
have a repetitive structure of arguments and presenta¬ 
tion, so that the reader can choose to read the first sub¬ 
section, skip ahead to the Conclusions, and return to read 
Subsections B and C at leisure. 


A. Catalytic Janus particles with high coverage 
and uniform surface mobility 

Previously m we studied the dynamics of a model 
catalytically active Janus particle suspended in a quies¬ 
cent fluid and near a wall. We found that, for a uniform 
surface mobility 6, in the course of time a particle with 
a very high coverage by catalyst can be stably attracted 
to a “hovering” state in which it remains motionless at a 
height h* and angle £* = 3n/2 (see Fig. |2jb)), i .e., with 
its catalytic cap oriented vertically and away from the 
wall ED- In this state, depicted in the left panel of Fig. 
[4j the tendency of the particle to translate away from its 
cap (i.e., trying to avoid high solute concentrations) is 
balanced by the accumulation of a solute “cushion” near 
the impenetrable wall (which is due to the confinement 
of the solute between the particle and the wall). The sta¬ 
bility of this hovering state against perturbations in the z 
direction can be understood easily: if the particle moves 
closer to the wall, solute accumulation is enhanced, and 
repulsion from the wall is stronger; if the particle moves 
away from the wall, repulsion from the wall is weakened, 
but the particle still translates away from its cap. Less 
obviously, hydrodynamic interactions with the wall stabi¬ 
lize the particle against perturbations in the angle £. Hy- 
drodynamically, the particle is characterized as a “puller” 
(see the flow lines in Fig. 3(b) in Ref. (HD)- Pullers are 
known to orient themselves perpendicular to planar sur¬ 
faces via hydrodynamic interactions pH] J^j 

The phase plane for a hoverer with xo = ~ cos i/j = 
0.85 and uniform surface mobility b in a quiescent fluid 


3 Additionally, while our investigation of hoverers initially assumed 
uniform surface mobility, we have found that these mechanisms 
are preserved if the surface mobilities on the cap and on the inert 
regions are unequal, i.e., for a wide range of bi ner t/b ca p A 1. 
Therefore the mechanism discussed here holds more generally. 
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(a) no external flow 



% / n 


(b) with shear flow 
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FIG. 5. (a) Phase plane for a Janus particle with high coverage by catalyst (xo = 0.85) and uniform surface mobility if the 

director lies in the yz plane ( p x = 0) and if there is no external flow. The filled magenta circle indicates an attractive “hovering” 
state. In this state, the particle remains motionless at a height h*/R = 1.21, with its cap oriented vertically and away from 
the wall (£* = 37 t/20 = 270°). The tendency of the particle to move away from its cap is balanced by the accumulation of 
solute near the wall. Clearly, the basin of attraction for this state is large. Additionally, many trajectories (roughly speaking, 
the whole basin of attraction) converge without overlapping to a “slow manifold” (the two visible branches merging into the 
attractor), indicative of two timescale dynamics: the vertical motion is much faster than the rotation of the particle [31] . 
Trajectories with an initial angle 0 < < 1 escape from the wall. The mirror symmetry of the region 1 < £/n < 2 across 

£/7t = 3/2, as well as the mirror symmetry of the region 0 < £/n < 1 across (;/n = 1/2, is due to rotational symmetry around 
the z axis in the absence of flow, (b) Phase plane for the same Janus particle if the director lies in the shear plane (p x = 0) 
and in the presence of a shear flow with yR/Uo = 0.02. Note that the phase plane is periodic in the £ direction. Due to the 
rotation of the particle by the flow, the attractor (filled green circle) has migrated into the region 1 < £/-7r < 3/2, and is now 
at h*/R = 1.35 and £* = 257°. As discussed in Sec. II, this is the region in which positive rheotaxis is possible. A particle in 
this state moves upstream with a steady height and orientation. The orange cross indicates a saddle point. An enlarged view 
of the region containing the saddle point and the attractor is provided in the right panel; some trajectories have been omitted 
for clarity. 


(i.e., no external flow) is shown in Fig. [5j[a) . The phase 
plane indicates the evolution of the particle configuration 
(h 1 £) for any initial condition (/io,£o)- I n the absence 
of flow, the system is symmetric for rotations around 
the z axis. This rotational symmetry is evinced by two 
mirror symmetries in the phase plane: symmetry of the 
region 0 < £/ 7 r < 1 for reflection across £/ 7 r = 1/2, 
and symmetry of the region 1 < £/tt < 2 for reflection 
across £/ir = 3/2. The “hovering” attractor is shown as 
a filled magenta circle at £* = 3tt/2 and h* /R = 1 . 21 . 
As discussed in Ref. [3T], many trajectories in the basin 
of attraction are drawn to a “slow manifold” indicative 
of two timescale dynamics: rotation is much slower than 
motion in z direction. For initial conditions £o < 7 r, the 
particle escapes the surface. 

Now we consider the same “hoverer” in shear flow. The 
flow will rotate the particle clockwise, decreasing £ from 
£ = 37 r/ 2 . On the other hand, the hydrodynamic in¬ 
teraction which stabilized “hovering” will tend to rotate 
the particle back towards £ = 37 r/ 2 . As depicted in the 
center and right panels of Fig. [4j these two contribu¬ 
tions to rotation can balance for sufficiently large angular 
displacement from the cap-down orientation. Hence the 
right panel of Fig. [4] provides a potential realization of 


the mechanism illustrated in Fig. [2] Therefore we an¬ 
ticipate that, at least for sufficiently slow external flows, 
there will a stable angle tv < < 37 t / 2 . Moreover, we 

anticipate that the “cushion” effect will be preserved to 
produce a steady height h*. The numerical simulations 
for, e.g., yR/Ug = 0.02 lead to a phase plane as shown in 
Fig. [5jb) and confirm these expectations: the attractor 
(filled green circle) is preserved and it migrates into the 
region 1 < £/ 7 r < 3/2, being now located at h*/R = 1.35 
and £* = 257°. Trajectories from a large section of phase 
space are drawn to this attractor. (Note that the phase 
plane is periodic in the £ direction.) Additionally, we 
indicate a saddle point (orange cross). The saddle point 
and attractor are rather close, as we have chosen a shear 
rate close to the (numerically estimated) upper critical 
value 7 C ss 0.021 x Uq/R. At the critical value, the sad¬ 
dle point and the attractor collide and annihilate each 
other. Above the critical shear rate, there is no stable 
rheotaxis, and, for this choice of the surface chemistry, 
all trajectories escape from the surface. We have chosen 
yR/Uo near the upper critical shear rate since, as will 
be discussed in the Conclusions, “strong” external flows 
have the greatest experimental relevance and accessibil¬ 
ity. Additionally, the relaxation time scale in Eq. (10) 
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FIG. 6. Trajectories and phase maps for a catalytic Janus particle with high coverage by catalyst (xo = 0.85) and uniform 
surface mobility. A trajectory starts from an initial height ho, an initial director orientation do and (j> o, and an initial lateral 
position ( xo,yo ) = (0,0). There is an ambient shear flow with dimensionless shear rate a /R/Uq = 0.02. The panels (a), (b), 
and (c) show trajectories launched from ho/R = 2, ho/R = 4, and ho/R = 6, respectively. Green trajectories are attracted to 
a rheotactic swimming state, whereas blue trajectories escape from the surface. In each panel, a white circle shows the initial 
spatial position of all trajectories. Arrowheads on selected trajectories indicate the direction of motion. In panels (d), (e), and 
(f) we show phase maps which indicate the behavior for each initial height and orientation. The filled green squares correspond 
to the rheotactic trajectories, and the open blue circles correspond to the escaping trajectories. Note that, for visual clarity, 
not all of the trajectories considered in (d) through (f) are plotted in (a) through (c). In particular, escaping trajectories are 
only shown in the region x > 0 . 


decreases as the shear rate is increased. Therefore, we ex¬ 
pect the approach towards the rheotactic state to occur 
most rapidly for 7 R/U 0 being close to the upper critical 
shear rate. Although we leave an exhaustive parametric 
study to future research, we note that results of addi¬ 
tional numerical calculations at lower values of 7 R/Uq , 
omitted here, agree with this expectation. 

Since the attractor is in the region tt < £* < 37 r/ 2 , 
it should be stable against perturbations of the director 
out of the shear plane, as discussed in Sec. II. In order to 
probe the stability of the attractor and its basin of attrac¬ 
tion, we consider an ensemble of trajectories launched 
from various initial director angles 9q and </>o (see Fig. 
[lja)) , initial heights ho/R = 2 , h 0 /R = 4, and ho/R = 6 , 
and the initial lateral position (a:o, yo ) = (0,0). For these 
three initial heights Figs. (b), and (c) show three- 

dimensional trajectories of the coordinates (x, y, z) of the 
centers of the particles. Each trajectory has been ob¬ 
tained from an initial position and orientation by numer¬ 
ically integrating Eqs. (U)-©. Green trajectories are 
“rheotactic.” Particles fo lowing these trajectories are 
attracted to the steady height h* /R = 1.35 and angle 


£* = 257° and move upstream. Particles following blue 
trajectories escape from the surface (i.e., the trajectories 
cross h/R = 7.1 from below.) Phase maps, indicating 
how the particle behavior depends on the initial height 
and orientation, are shown in Figs. id), (e), and (f), 
which correspond to ho/R = 2 , h 0 /R = 4, and ho/R = 6 , 
respectively. Clearly, rheotaxis is achieved for a large 
basin of initial conditions. 

An example of a rheotactic trajectory is shown in 
Fig. 0 Starting from the initial orientation ho/R = 6 , 
60 = 120°, and </>o = 67.5°, the particle has nearly at¬ 
tained the steady height and orientation after moving 
only a few particle diameters. The particle is attracted 
to a configuration in which it is tilted slightly away from 
the “hovering” state by the shear flow. As the (blue) cap 
is oriented slightly downstream (i.e., £* = 257° < 270°), 
the particle moves upstream. In Fig. [8j the particle 
starts near the wall, but pointing away from it. Due to 
this initial orientation, the particle moves a few diame¬ 
ters away from the wall, where the hydrodynamic and 
chemical influence of the wall is very weak. However, the 
flow rotates the particle to swim back towards the wall. 
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FIG. 7. Trajectory for a catalytic Janus particle with high coverage by catalyst (xo = 0.85, blue) and uniform surface mobility 
for the initial conditions ho/R = 6 , do = 120°, and <j >o = 67.5° in a shear flow with 7 R/Uo = 0.02. The particle is attracted to a 
configuration in which it is tilted slightly away from the “hovering” state (i.e., a completely vertical orientation) by the flow, and 
consequently moves upstream. Particle positions along the trajectory are labeled by the dimensionless times t/To at which they 
occur. The particle is clearly almost in the rheotactic state at t/To = 250, which we estimate in the Conclusions to correspond 
to t ~ 30 s for typical catalytic Janus particles used in experiments. Note that we have slightly enlarged the appearance of the 
inert part of the particle surface for visual clarity. For this trajectory, the time dependence of various quantities is shown by 
solid lines with stars in Fig. [9] 



FIG. 8 . Trajectory for a catalytic Janus particle with high coverage by catalyst (xo = 0.85, blue) and uniform surface mobility 
for the initial conditions ho/R = 2 , 6q = 40°, and (po = 135° in a shear flow with r yR/Uo = 0 . 02 . The particle initially 
moves away from the wall, owing to the orientation of its catalytic cap. A few particle diameters away from the wall, the 
influence of the wall is very weak. However, the particle is rotated by the flow to move back towards the wall. As the particle 
approaches the wall, the influence of the wall strengthens. The particle rotates into the “tilted hoverer” configuration and 
moves upstream. Particle positions along the trajectory are labeled by the dimensionless times t/To at which they occur. Note 
that we have slightly enlarged the appearance of the inert part of the particle surface for visual clarity. For this trajectory, the 
time dependence of various quantities is shown by solid lines with squares in Fig. [9] 


Upon returning to the vicinity of the wall, the particle 
rotates into the “tilted hoverer” configuration. 

For selected rheotactic trajectories, in Fig. [9] we plot 


the time evolution of the particle height and of the com¬ 
ponents of the director. As expected, the component p x 
of the director perpendicular to the shear plane decays 
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FIG. 9. (a)-(d) Time evolution of various quantities for five of the rheotactic trajectories of the catalytic Janus particle 

considered in Fig. [6] The quantity p x decays to zero, while p y , p z , and h attain asymptotically non-zero values p* = —0.22, 
pi = —0.98, and n/R = 1.35. The initial conditions for these trajectories are: solid line with stars, ho/R = 6, do = 120°, 
(j> 0 = 67.5° (see also Fig. [7j; solid line with filled circles, ho/R = 2, do = 180°, 4>o = 0°; dotted line, ho/R = 4, do = 130°, 
4>o = 135°; dashed line, ho/R = 6, do = 70°, tj >o = 112.5°; solid line with squares, ho/R = 2, do = 40°, <f>o = 135° (see also Fig. 
[8|. In panel (e), \p x \ is plotted for all studied rheotactic trajectories on a semi-logarithmic scale. The exponential character of 
the decay of \p x \ is clearly visible. The slope of the dotted black line indicates the exponential time dependence predicted by 
Eq. :fm . 


to zero. The height h and the two in-plane components 

t and p z attain asymptotically non-zero values. In Fig. 

e), we plot the time evolution of \p x \ for all rheotac¬ 
tic trajectories studied for Fig. [bj The decay of \p x \ is 
clearly exponential, and the timescale for decay closely 
agrees with the prediction of Eq. ( flO] ) , which is plotted 
as the dotted line in Fig. [9je) . Note that in Figs. [7J |8j 
and [9] time is given in dimensionless units as t/To- In a 
previous section, T 0 has been defined as T 0 = Uq/R. As 
will be discussed in the Conclusions, we estimate T 0 to 
be T 0 « 0.125 s for catalytic Janus particles as typically 
used in experiments. 


B. Catalytic Janus particle: half coverage and 
inhomogeneous surface mobilities 

In order to demonstrate the general character of our 
theoretical findings, we seek to reach the rheotactic state 
of Fig. §b) via the alternative pathway of a different 
surface chemistry designed to realize a distinct physical 
mechanism. Specifically, we consider a particle which is 
half covered by catalyst (xo = 0 ), but the surface mo¬ 
bilities of the inert region and of the catalytic cap are 
taken to differ: bi ner t 7 ^ b cap . In part, this choice is mo¬ 
tivated by the fact that by default experimental studies 
use particles with half coverage, because these can read¬ 


ily be prepared by vapor deposition. Moreover, since 
the catalytic and inert surface regions consist of different 
materials, they are likely to give rise to different surface 
mobilities, too. 

In our previous study of systems without flow, we have 
isolated two distinct wall-induced contributions to the ro¬ 
tation of a Janus particle pITi . One contribution is due 
to the hydrodynamic interaction of the particle with the 
wall. Disturbance flows created by the motion of the par¬ 
ticle are reflected by the wall, coupling back to the par¬ 
ticle. We found that, for half coverage (xo = 0), hydro- 
dynamic interactions always tend to rotate the catalytic 
cap of the particle towards the wall (and thus the direc¬ 
tor away from the wall). Therefore, this contribution to 
particle rotation cannot oppose the rotation by the shear 
flow (cyan arrow in Fig. |2jb)) for 7 T < £ < 37 r/ 2 , i.e., it 
cannot provide the magenta arrow in Fig. Sb). How¬ 
ever, we also found that if b inert 7 ^ b capi wall-induced 
chemical gradients can contribute to particle rotation. If 
binert/b C ap < 1 , repulsion of the solute (i.e., the reac¬ 
tion product) from the inert region is weaker than repul¬ 
sion from the catalytic cap. Accordingly, wall-induced 
chemical gradients (i.e., a higher concentration of solutes 
on the side of the particle surface closer to the wall) 
tend to rotate the catalytically active cap away from the 
wall. (Note that chemical gradients will not drive rota¬ 
tion of the particle in the absence of the wall because the 






















14 


(a) no external flow 



(b) with shear flow 



FIG. 10. (a) Close-up of the phase plane of a Janus particle with half coverage (\o = 0) and inhomogeneous surface mobilities 

binert/bcap = 0.3 if the director is oriented in the yz plane and there is no external flow. There is an attractor (filled magenta 
circle) at h*/R = 1.03 and £* = 204°. At this point, the hydrodynamic interaction with the wall, which tends to rotate the 
cap towards the wall, is balanced by the effect of wall-induced chemical gradients, which tend to rotate the cap away from the 
wall. Importantly for the behavior in the presence of flow (see (b)), there is a region in which the net effect of the activity 
of the particle is to rotate the cap away from the wall. Within the region 1 < £/-7r < 3/2, this occurs where the phase space 
trajectories in (a) flow towards larger £ (i.e., to the right.) (b) Phase plane for dynamics in the shear plane ( p x = 0) of the 
same particle with an external shear flow r yR/Uo = 0.05. There is a rheotactic attractor (filled green circle) at h*/ R = 1.08 
and = 200°, as well as a saddle point (orange cross). At the attractor, the net orientational effect of the chemical activity, 
which tends to rotate the cap away from the wall, is balanced by the orientational effect of the shear, which tends to rotate the 
cap towards the wall, as illustrated in Fig. §b). Note that the phase plane is periodic in £. An enlarged view of the region 
containing the saddle point and the attractor is provided in the right panel. 


Janus particle is axially symmetric.) Therefore, taking 
binert/bcap < 1 generates a contribution to the rotational 
velocity of the particle which corresponds to the magenta 
arrow in Fig. UN- 

Therefore, we consider a particle with half coverage 
(Xo = 0) and bi nert /b cap = 0.3. A close-up of the phase 
plane for the dynamics is shown in Fig. for the 

case that the director lies in the yz plane and that there 
is no external flow. There is an attractor very close to 
the wall at h*/R = 1.03 and £* = 204°. At this point, 
the hydrodynamic interaction with the wall, which tends 
to rotate the cap of the particle towards the wall, is bal¬ 
anced by the effect generated by wall-induced chemical 
gradients, which tends to rotate the cap away from the 
wall. The particle moves along the wall with a steady 
height and steady orientation, which we called a “sliding” 
steady state (311 ■ Importantly, there are regions of the 
phase space in which the rotational effect of wall-induced 
chemical gradients is stronger than the rotational effect 
of the hydrodynamic interaction, such that in sum the 
cap tends to rotate away from the wall. In particular, in 
the interval 1 < //n < 3/2 in Fig. HU*), the cap rotates 
away from the wall whenever trajectories flow towards 
larger £ (see Fig. |2j(b)), i.e., to the right of the plot in 
Fig. |10| (a). In this region, rotation away from the wall 
by chemical activity is so strong that it can balance the 
rotational effect of shear towards the wall, as illustrated 
in Fig. [2](b). In addition, h = 0 whenever the tangent 


to the trajectories h(£) in Fig. [lOja) is horizontal, i.e., 
d[h(£(t))\/dt = £ ( dh/d £) = 0. Therefore, the swimming 
activity of the particle can potentially also on its own 
produce a steady height h *. We also note that some tra¬ 
jectories cross below the minimum height h/R = 1.02 
which we consider. The particles seemingly “crash” into 
the wall. It is likely that many of these trajectories re¬ 
cross the line h/R = 1.02 and flow to the attractor if the 
numerical calculations were extended below h/R = 1.02; 
however, as discussed above, the physical approximations 
inherent in the calculations break down very close to the 
wall0 

We now consider the same particle in shear flow. The 


4 Without any major difficulties, the numerical calculations, 
within the current mathematical description, can be extended 
technically to separations below h/R < 1.02, corresponding to 
a particle/wall gap of less than 50 nm for a particle of radius 
R = 2.5/mi. However, these results would be physically ques¬ 
tionable because our model cannot be expected to apply in that 
range. As discussed previously, for small particle/wall separa¬ 
tions one can no longer assume a separation of length scales be¬ 
tween solute/particle interactions and bulk hydrodynamic flow. 
Moreover, in that range, the typical surface interactions (such 
as van der Waals and electrostatic double layer forces) are rel¬ 
evant and therefore the assumption that the active particle is 
force and torque free breaks down. Such contributions can be 
included within a more complex model of active particles, which 
can be studied along similar lines. 
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FIG. 11. Trajectories and phase maps for a catalytic Janus particle with half coverage (xo = 0) and inhomogeneous surface 
mobilities bi n£r t/b ca p = 0.3. A trajectory starts from an initial height ho, an initial director orientation do and cj> o, and the 
initial lateral position ( xo,yo) = (0,0). There is an external shear flow with dimensionless shear rate 7 R/Uo = 0.05. The 
panels (a), (b), and (c) show trajectories launched from ho/R — 2, ho/R = 4, and ho/R = 6 , respectively. Green trajectories 
are attracted to a rheotactic swimming state, blue trajectories escape from the surface, and red trajectories “crash” into the 
wall (i.e., they approach the surface closer than the minimal height we consider numerically). In each panel, a white circle 
shows the initial spatial position of all trajectories. Arrowheads on selected trajectories indicate the direction of motion. In 
panels (d), (e), and (f) we show phase maps which indicate the behavior for each initial height and orientation. The filled 
green squares correspond to the rheotactic trajectories; the open blue circles correspond to the escaping trajectories; and the 
red crosses correspond to the “crashing” trajectories. Note that, for visual clarity, not all of the trajectories considered in (d) 
through (f) are plotted in (a) through (c). 


phase plane with shear rate jR/Uo = 0.05 is shown in 
Fig. [Tojb). As anticipated, there is a rheotactic attractor 
(filled green circle), which is located at h*/R = 1.08 and 
£* = 200°. Additionally, there is a saddle point (orange 
cross). As in the case of the “hoverer,” we have chosen a 
shear rate close to the upper critical value (numerically 
estimated to be , y c R/Uo ~ 0.06 for this surface chem¬ 
istry), and therefore the saddle point and attractor are 
very close to each other. 

As for Fig. [6j we launch an ensemble of particles from 
the lateral position ( Xo,yo ) = (0,0) of the center of the 
particles, for various initial director orientations 9 0 and 
(f> 0 , and heights ho/R = 2, ho/R = 4, and ho/R = 6. 
The three-dimensional trajectories of the center of the 
particles are shown in Figs. H3a), (b), and (c). The 
green trajectories are rheotactic. Particles following blue 
trajectories escape from the surface. Finally, the red tra¬ 
jectories are those which “crash” into the surface, as dis¬ 
cussed above. Phase maps indicating the particle behav¬ 
ior as a function of the initial orientation are shown in 
Figs. ©d), (e), and (f). For selected rheotactic trajec¬ 


tories, in Fig. [l2|(a)-(d) we plot the time evolution of the 
height h and of the director components p x , p y , and p z . 
In Fig. 12(e), we plot \p x \ for all rheotactic trajectories 
studied in Fig. 11 As in Fig. [9je) , the decay time for 
p x closely agrees with the prediction of Eq. (101, plotted 
as the dotted line in Fig. 12 (e). Finally, a representative 
rheotactic trajectory is shown in detail in Fig. 13 The 
radius of curvature of the trajectory in its evolution to¬ 
wards the rheotactic steady state is clearly much larger 
than for the “hoverer” shown in Fig. [7] As with the 


hoverer, the dimensionless times t/To in Figs. 12 and 13 


can be converted to dimensional, and thus experimentally 
relevant values via the estimate T 0 « 0.125 s. 


C. Squirmer 

In order to further reveal the general character of our 
theoretical results and predictions, we consider a spheri¬ 
cal “squirmer.” A squirmer has a prescribed surface ve¬ 
locity. It interacts with a bounding wall hydrodynami- 
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t/T 0 


FIG. 12. (a)-(d) Time evolution of p x , p y , p z , and h/R for five of the rheotactic trajectories considered in Fig. 11 The 

quantity p x decays to zero, while p y , p z , and h attain asymptotically non-zero values p* = —0.94, pi = —0.35, and h*/R = 1.08. 
The initial conditions for these trajectories are: solid line with stars, ho/R = 6, do = 160°, (f> o = 225°; solid line with filled 
circles, ho/R = 2 , do = 130°, <j> o = 315°; dotted line, ho/R = 4, do = 60°, (j>o = 0°; dashed line, ho/R = 6, d 0 = 90°, <j> o = 22.5° 
(see also Fig. 131; solid line with squares, ho/R = 2, do = 70°, 4 >o = 180°. In panel (e), \p x \ is plotted for all studied rheotactic 
trajectories on a semi-logarithmic scale. The asymptotic decay of \p x 
indicates the exponential time dependence predicted by Eq. |To|. 


is clearly exponential. The slope of the dotted line 



FIG. 13. Trajectory for a catalytic Janus particle with half coverage (xo = 0) and inhomogeneous surface mobilities 
binert/b C ap = 0.3 for the initial conditions ho/R = 6, do = 90°, and </>o = 22.5° in a shear flow with jR/Uo = 0.05. Par¬ 
ticle positions along the trajectory are labeled by the dimensionless times t/To at which they occur. The particle has “turned 
around” at t/To = 350, which we estimate in the Conclusions to correspond to t ~ 45 s for typical catalytic Janus particles 
used in experiments. For this trajectory, the time dependence of various quantities is shown by dashed lines in Fig. |12| 


cally, but not chemically, because the surface flow is not 
driven by a distribution of solute. The prescribed sur¬ 
face velocity is often taken to model the time-averaged 
motion of cilia on the surface of a micro-organism. 

We follow Li and Ardekani [35] in restricting our con¬ 
sideration to the first two squirming modes, given by B i 


and B 2 in Eq. (11). Li and Ardekani found that stable 
dynamical attractors exist for various values of B 2 /B 1 
and of the Reynolds number Re [38]. (The lowest value 
of Re they considered was Re = 0.1.) For instance, 
for B 2 /B 1 = 3 and Re = 1, they found an attractor 
with h*/R = 1.47 and 6 * = 103.2°. By construction 
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FIG. 14. (a) Phase plane for a squirmer with B 2 /B 1 = 7 (Eq. 0) if the particle director lies in the yz plane and if there is 

no external flow. There are stable dynamical attractors (filled magenta circles) for which the particle achieves a steady height 
and a steady orientation and “slides” along the wall. For these attractors, the particle director points towards the wall. There 
are also unstable fixed points (open magenta circles) and one saddle point (orange cross). Due to rotational symmetry around 
the z direction and in the absence of flow, the region 0 < £/n < 1 is mirror symmetric across the line = 1/2, and the region 
1 < C/ 71 " < 2 is mirror symmetric across the line £/n = 3/2. (b) Phase plane for the same squirmer if the particle director lies 
in the yz plane and if there is a shear flow of strength 'yR/Bi = 0.09. There is a rheotactic attractor (filled green circle). Since 
it lies in the region 1 < £/n < 3/2, this attractor is stable against small perturbations of the director out of the yz plane (i.e., 
for p x ^ 0), as predicted in Sec. II. There is another fixed point (filled magenta circle) which is stable against perturbations 
within the plane p x = 0, but unstable against perturbations out of this plane. Likewise, this instability has been predicted in 
Sec. II for fixed points in the region 3/2 < £/7r < 2. There are two saddle points (orange crosses). The phase plane is periodic 
in the £ direction. An enlarged view of the region containing the rheotactic attractor and the nearby saddle point is provided 
in the right panel. 


our numerical method probes the limit Re = 0. For 
B 2 /B 1 = 3 (and Re = 0), we have found an attractor 
at h*/R = 1.64 and 6 * = 102.8°. We note that the 
dynamics of a squirmer near a boundary, including the 
possibility of moving at a stable height and orientation, 
was studied also by Ishimoto and Gaffney for Re = 0 E3- 
Here, we consider a squirmer with B 2 /Bi = 7. We 
choose this large ratio in order to achieve a strong hy¬ 
drodynamic interaction with the wall, permitting rheo- 
taxis for a relatively high shear rate "fR/Bi = 0.09, as 
will be discussed below. (By comparison, for B 2 /B\ =3 
rheotactic states occur only for jR/Bi < 0.002.) In Fig. 
14 'a) we show the phase plane for B 2 /Bi = 7 if the 
director lies in the yz plane and if there is no external 
flow. There are two stable dynamical attractors (filled 
magenta circles) for which the particle moves at a steady 
height and a steady orientation with its director point¬ 
ing towards the wall, similar to the “sliding” states we 
found for catalytic Janus particles pH] . Due to the ro¬ 
tational symmetry of subproblem II (see Fig. [ 5 ] and the 
corresponding discussion in Sec. III.A), these two attrac¬ 
tors are actually the same in the sense that they corre¬ 
spond to the same sliding state, only that the particle 
slides towards the negative and positive y direction for 
the attractor with £*/7 t < 3/2 and £*/7r > 3/2, respec¬ 
tively. In addition, there are two unstable fixed points 
with £*/7 t < 1 (open magenta circles; this is again the 


same fixed point by symmetry) and a saddle point (or¬ 
ange cross). 

Now we consider the effect of an external flow with 
'yR/Bi = 0.09 on the two attractors. They remain 
attractors for in-plane dynamics and stay at approxi¬ 
mately the same locations in the phase plane, as can be 
seen in Fig. 0b). The fixed point at h*/R = 1.22 
and £* = 207° (filled green circle) lies in the region 
1/2 < £ /7T < 3/2 and therefore within the interval of 
£ for stable rheotaxis. In this configuration, the direc¬ 
tor is oriented upstream and towards the wall, as in Fig. 
§b). This fixed point is therefore a global attractor or an 
attractor for fully three-dimensional dynamics. For the 
other fixed point at h*/R = 1.14 and £ = 332° (filled ma¬ 
genta circle), the particle director is oriented downstream 
(p v > 0) and towards the wall (p z < 0); according to Eq. 
(101, this fixed point should be linearly unstable against 


perturbations out of the plane, i.e., away from p x = 0. 
In addition, there are two saddle points (orange crosses). 

As for Figs. [6] and 11 we launch an ensemble of parti¬ 
cles from the lateral position (xo, yo) = (0,0) of the center 
of the particles, for various initial angles 9q and <f>o, and 
heights ho/R = 2, h 0 /R = 4, and h 0 /R = 6. The three- 
dimensional trajectories of these particles are shown in 
Figs. 0a), (b), and (c). Green indicates rheotactic tra¬ 
jectories; blue indicates trajectories which escape from 
the surface; and red trajectories “crash” into the wall 
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FIG. 15. Trajectories and phase maps for a squirmer with B2/B1 = 7 (Eq. (111). A trajectory starts from an initial 
height ho, an initial director orientation 9 o and <j> 0 , and the initial lateral position ( xo,yo) = (0,0). There is an external shear 
flow with dimensionless shear rate -yR/Bi = 0.09. The panels (a), (b), and (c) show trajectories launched from ho/R = 2, 
ho/R = 4, and ho/R = 6, respectively. Green trajectories are attracted to a rheotactic swimming state; blue trajectories 
escape from the surface; and red trajectories “crash” into the wall (i.e., they approach the surface closer than the minimal 
height we consider numerically). In each panel, a white circle shows the initial spatial position of all trajectories. Arrowheads 
on selected trajectories indicate the direction of motion. In panels (d), (e), and (f) we show phase maps which indicate the 
behavior for each initial height and orientation. The filled green squares correspond to rheotactic trajectories; the open blue 
circles correspond to the escaping trajectories; and the red crosses correspond to the “crashing” trajectories. Note that, for 
visual clarity, not all of the trajectories considered in (d) through (f) are plotted in (a) through (c). 


(i.e., they approach the wall closer than h/R = 1.02, 
which is the minimum height we consider numerically). 
The phase maps in Figs. [Hd), (e), and (f) indicate the 
particle behavior as a function of the initial height and 
the initial orientation. One rheotactic trajectory is shown 


in detail in Fig. 16 


For selected rheotactic trajectories, in Figs. Hla)-(d) 
we plot the time evolution of the height and of the ori¬ 
entation of the particle. Since the rheotactic attractor 
is oscillatory (see Fig. jdjb)), these quantities exhibit 
damped oscillations. The time evolution of \p x \ for all 
rheotactic trajectories studied is plotted in Fig. &). 
As for the catalytically active particles, the exponential 
decay of the component p x closely agrees with the one 
predicted by Eq. (10); in Fig. [T7|(e) this theoretical 


prediction is indicated by the slope of the dotted black 
line. For the squirmer, we can test another prediction 
of the linear stability analysis. As discussed previously 
in Fig. H3 b ), there is a fixed point, represented by the 
filled magenta circle, which is, as obtained numerically, 
a stable attractor for in-plane dynamics, i.e., for p x = 0. 


At this fixed point, the particle moves downstream at a 
steady height and a steady orientation. However, this 
fixed point lies in a region (3/2 < £/7 t < 2) for which 
linear stability analysis predicts fixed points to be un¬ 
stable against small perturbations in p x . In this con¬ 
text we launch two trajectories with initial conditions 
near this fixed point (Fig. 


18). 


Both trajectories start 
with ho/R = 1.2 and 9q = 150°, but one with 4>q = 90° 
(dashed magenta line in Fig. |l8|(a)), i.e., p x = 0, and 
the other with </> 0 = 89° (solid green line), i.e., p x ^ 0. 
For the trajectory with <p 0 = 90°, the particle always 
remains in the shear plane ( p x = 0) and is attracted 
to the downstream-moving state. However, for the tra¬ 
jectory with </> 0 = 89°, the particle moves out of the 
shear plane, and is ultimately attracted to the rheotac¬ 
tic upstream-moving configuration. In Fig. HIb), we 
show the time evolution of \p x \ for the rheotactic tra¬ 
jectory. Both the initial instability and the asymptotic 
approach to the rheotactic state clearly behave exponen¬ 
tially. The time dependences predicted by Eq. ( flO| ) are 
illustrated by dotted black lines. We note that the time 
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FIG. 16. Rheotactic trajectory for a squirmer with B 2 /B 1 = 7 and initial condition ho/R = 6, do = 130°, and cj >0 = 225° in 
a shear flow with A /R/Bi = 0.09. The black shaded part of the particle indicates its rear. The trajectory exhibits oscillations. 
Selected particle positions along the trajectory are labeled by the dimensionless times t/T a at which they occur. For this 
trajectory, the time dependence of various quantities is shown by solid lines with stars in Fig. [IT] 



FIG. 17. (a)-(d) Time evolution (in units of T a = R/B 1 ) of various quantities for selected rheotactic trajectories of the 

squirmer discussed in Fig. [15] The quantity p x decays to zero, while p y , p z , and h attain asymptotically non-zero values 
Py = —0.89, pt = —0.45, and h*/R = 1.22. The initial conditions for the trajectories are: solid line with stars, ho/R = 6, 
6 <o = 130°, (f >0 = 225° (see also Fig. |16[ ); solid line with filled circles, ho/R = 4, do = 110°, <f >0 = 180°; dotted line, ho/R = 4, 
do = 120°, 4 >0 = 315°; dashed line, ho/R = 6, do = 110°, (j >0 = 0°; solid line with squares, ho/R = 2, do = 110°, <j)o = 157.5°. 
Note that, for visual clarity, in (c) only three trajectories are plotted. In (e), \p x \ is plotted on a semi-logarithmic scale for all 
rheotactic trajectories studied. The slope of the dotted line indicates the exponential time dependence predicted by Eq. (101. 


scales for initial growth and asymptotic decay of \p x | are 
similar. This can be inferred from Fig. M h ) , in which 
the filled magenta circle and the filled green circle rep¬ 
resent the downstream-moving and rheotactic upstream- 
moving fixed points, respectively. The original mirror 


symmetry across £/-7r = 3/2 of Fig. 14(a), in which there 
is no external flow, is approximately preserved, and hence 
the two fixed points have approximately the same values 
of |p* I, |p*|, and h*. (As a reminder, Fig. J2j(a) relates p 
and f.) Since the time scale for growth or decay of A p x is 
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determined by |p*|, |p*|, and the local flow rate 7 f(h*/R) 
(see Eq. 10), the two time scales are approximately the 
same for the two fixed points. 


IV. CONCLUSIONS 

We have theoretically investigated the possibility that 
a spherical active particle with propulsion mechanism, 
which (i) is axisymmetric and (ii) can be described in 
terms of an effective slip velocity, may exhibit rheotaxis 
(and, in particular, upstream rheotaxis) in shear flow 
near a planar surface. (We define rheotaxis as to de¬ 
note the approach of the particle to a robust and stable 
steady state in which the orientation vector lies within 
the plane of shear). We have found that rheotaxis of such 
a particle is indeed possible, even though its spherical 
geometry rules out the intuitive “weather vane” mecha¬ 
nism generally invoked in order to explain rheotaxis for 
elongated microswimmers (e.g., sperm, E. coli , and poly¬ 
mer/hematite dimers.) 

Furthermore, we have shown that for any such parti¬ 
cles it is sufficient to analyze the dynamics of the particle 
with the orientation vector lying in the shear plane in or¬ 
der to determine whether or not rheotaxis occurs. In 
particular, for positive (upstream) rheotaxis, there must 
be an attractive steady state for the in-plane dynamics 
in which the orientation vector points upstream and to¬ 
wards the surface. For such a configuration, any small 
perturbation of the orientation vector out of the shear 
plane is damped by the moving activity of the particle. 
These findings significantly simplify the theoretical calcu¬ 
lations. Moreover, they provide conceptually simple rules 
for the design of artificial rheotactic microswimmers. In 
order to achieve rheotaxis, one only needs to tailor the 
self-propulsion mechanism of the particle such that two 
criteria are satisfied: (i) the contribution of near-surface 
motion to the rotation of the particle can stably bal¬ 
ance shear-induced rotation, as shown in Fig. §b), and 
(ii) the particle, through its near-surface moving activ¬ 
ity, can achieve a steady state of fixed height. Specific 
examples of design for rheotaxis have been provided by 
exploring two distinct pathways for fulfilling these crite¬ 
ria via tailoring the surface chemistry of a catalytically 
active Janus particle. In the first pathway, the phoretic 
mobility is distributed homogeneously over the surface 
of the active particle, and the coverage by catalyst is ad¬ 
justed to provide the required activity-induced rotation. 
In the second pathway, the coverage by catalyst is kept 
to the value of one half, which is simpler for experimental 
realizations, but the phoretic mobility is taken to be dis¬ 
tributed inhomogeneously over the surface of the active 
particle. In both these cases, the numerical solutions of 
the equations of motion evidence rheotaxis and confirm 
the analytical prediction for the decay time for the com¬ 
ponent of the orientation vector out of the shear plane. 

The numerical study we presented in Sec. |III| A and B 
revealed the existence of a threshold value r c of the di¬ 


mensionless parameter r = 7 R/Uq above which rheotaxis 
is no longer possible. It is therefore important to under¬ 
stand the consequences of this constraint in the context 
of experimental studies with catalytically active Janus 
particles. We can judge the feasibility of experimentally 
realizing rheotactic Janus active particles by estimating 
the values of the shear rates corresponding to the thresh¬ 
old values r c = 0.021 and r c = 0.06 for “hoverer” (Sec. 
|III| A) and “slider” (Sec. |III| B), respectively, which have 
been determined by our numerical calculations. Consid¬ 
ering, as in Sec. 0 E, a typical Janus particle of radius 
i? « 2.5 pm which is half-covered with catalyst and in 
the free space moves with a velocity U/. s . « 5 pm/s, 
i.e., Uq = 4 x Uf . s . ~ 20 pm/s 341] , we find that the 
maximal shear rate, above which rheotaxis is no longer 
possible, is 7 C = r c Uo/R = 0.17 s -1 for a “hoverer” and 
7 C = 0.48 s ” 1 for a “slider”, respectively. These values 
are low, but within or nearly within the range 7 = 0.2 s _1 
to 7 = 9 s _1 explored in Ref. [ 8 ], and within the range 
7 = 0.1 s -1 to 7 = 20 s -1 explored in Ref. f20j. In order 
to relax these bounds on the shear rate, the ratio R/Uq 
must be increased. We suggest two possible approaches 
to achieve this. The first one is to exploit the nonlin¬ 
ear relationship between size and propulsion velocity in 
the range of small particle radii, which was reported by 
Ebbens et al. @3]. Following Ref. [13), a 1 pm ra¬ 
dius particle with half coverage has a moving velocity of 
Uf . s . = 9 pm/s 031. For such an active Janus particle, 
half-covered by catalyst, the maximal shear rate, which 
corresponds to a slider and allows for rheotaxis, becomes 
7 C ~ 2.2 s -1 . However, this estimate should be consid¬ 
ered with due care. As the size of the active particle de¬ 
creases, the effects of the thermal noise on the orientation 
of the particle, which have been neglected in the present 
study, become significant and therefore the assumptions 
of the model leading to the predictions of the values for 
r c will break down. A second approach would be to ex¬ 
ploit the fact that the velocity scale Uq increases with the 
fuel concentration (at least within a certain range); this 
suggests using a high fuel concentration. For instance, 
Baraban et al. report that a R = 2.5 pm active Janus 
particle with half coverage can achieve speeds Uf . s . of 
more than 8 pm/s at a high concentration (15% vol¬ 
ume fraction) of hydrogen peroxide 04] . In this case, the 
maximal shear rate, which corresponds to a slider and 
at which rheotaxis can occur, would be in the range of 
0.8 s- 1 . 

Having estimated dimensional values of Uq and R for a 
catalytic particle, we obtain the time scale T 0 = R/Uq = 
0.125 s. We can acquire a rough sense of the time re¬ 
quired for rheotaxis in an experiment by converting the 
dimensionless times t/Tg given by the black numbers in 
Figs. 00 and [Ts] into dimensional quantities. For in¬ 
stance, for the “hoverer” in Fig. [7] the particle is clearly 
close to the rheotactic state by T/T 0 = 250. This cor¬ 
responds to an experimental time of ca 31s. Likewise, 
for the half-covered particle in Fig. [13] the particle has 
turned upstream by T/Tq « 350, corresponding to an 
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FIG. 18. (a ) Tw o trajectories with initial height and initial orientation near the fixed point represented by the filled magenta 

circle in Fig. |14[ b). This fixed point is a stable attractor for in-plane dynamics, i.e., for p x = 0, but it is unstable against small 
perturbations in p x . In this configuration, the particle moves in the downstream direction with a steady height and orientation. 
The two trajectories both have as initial conditions ho/R = 1.2 and do = 150°, as well as the same initial lateral position 
(*o, J/o) = (0, 0) (white circle). However, one trajectory (magenta dashed line) has <j> o = 90°, i.e., initially p x = 0, and therefore 
it is attracted to the downstream-moving state (i.e., the filled magenta circle in Fig. 14 b)). The other trajectory (solid green 
line) starts with tj> o = 89°, i.e., initially p x ^ 0, and ultimately it is attracted to the rheotactic upstream-moving state (i.e., the 
filled green circle in Fig. [Em (b) The time evolution of \p x \ in units of T a = R/B± for the rheotactic ( green) trajectory of 
(a), plotted on a semi-logarithmic scale. The slopes of the dotted lines indicate the exponential time dependences predicted by 
Eq. pi. 


experimental time of ca 44 s. Furthermore, for the two 
particle designs we can evaluate dimensional values of 
the relaxation time scale in Eq. (10 1 , which we denote 
as r a . For the hoverer, the dimensionless relaxation time 
is r a /T 0 = 25.4 (dotted line in Fig. |9je)), which ren¬ 
ders r a ss 3.2 s. Likewise, for the half-covered parti¬ 
cle, the dimensionless relaxation time is r a /To = 149 
(dotted line in Fig. [T2]e)), which renders r a « 19s. 
We can compare these values for r a with the time scale 
r r for rotational diffusion in order to estimate how ro¬ 
bust rheotaxis is against thermal noise. For a particle 
with R = 2.5 /Ltm in water at room temperature, the 
time scale r r for reorientation by rotational diffusion is 
DjT 1 = 87 j[iR?IksT ~ 95s [31]. Therefore, the rheotac¬ 
tic state is expected to be robust against thermal noise 
induced rotation out of the plane p x = 0 . 


Our theoretical findings are not restricted to catalyt- 
ically active Janus particles. In order to highlight their 
wide range of applicability, we have also investigated 
rheotaxis for a spherical “squirmer.” The squirmer model 
captures essential features of the motion of ciliated micro¬ 
organisms [45] and self-propelled liquid droplets [46]. Re¬ 
cently, it was shown theoretically that, for certain val¬ 
ues of the leading order squirming mode amplitudes, a 
squirmer can be attracted to a planar surface and can 
move at a steady height and orientation [37, 38]. In the 
presence of shear this attraction can allow one to fulfill 
the two criteria we have established for the occurrence 
of rheotaxis, and we have shown that rheotaxis does in¬ 
deed emerge. It is possible that the mechanism for rheo¬ 
taxis outlined in the present study is relevant for spher¬ 
ical or quasi-spherical micro-organisms, such as Volvox 


carteri. More speculatively, perhaps such organisms can 
dynamically adjust their squirming mode amplitudes to 
“turn on” and “turn off” rheotaxis. However, given the 
complexity of these organisms - Volvox , for instance, is 
bottom-heavy and spins around a fixed body axis |45] 
- experiments and more detailed numerical studies are 
needed to assess in this context the relevance of the mech¬ 
anisms studied here. Finally, we note that, while we focus 
on active particles which can be modeled via an effective 
slip velocity, we expect that our theoretical findings will 
be applicable to other spherical active particles with an 
axisymmetric propulsion mechanism, provided one can 
justify the linear separation into the two subproblems 
considered here of an active sphere in a quiescent fluid 
and an inert sphere in flow. 

A natural extension of our work would be to study 
elongated (e.g., ellipsoidal) active particles with an ax¬ 
isymmetric propulsion mechanism. The classic work of 
Jeffery demonstrates that inert ellipsoidal particles in un¬ 
bounded shear flow rotate in periodic orbits m- More 
recently, it was shown experimentally |48| and numeri¬ 
cally [49] that these orbits are preserved (although quan¬ 
titatively modified) if the particle is near, but not in steric 
contact with, a planar bounding surface. Theoretical and 
numerical studies of elongated active particles could shed 
light on how this orbital motion is transformed into the 
“weather vane” mechanism of rheotaxis with the addi¬ 
tion of active motion. An interesting point of compar¬ 
ison would be with Ref. [ 8 ], which presents a mathe¬ 
matical model for the rheotaxis of a spermatozoon in 
steric contact with a surface. Additionally, we anticipate 
that ellipsoidal artificial microswimmers would be able 
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to rheotax at higher shear rates than spherical particles, 
owing to the “weather vane” mechanism. Testing this 
expectation and quantifying the difference with spherical 
particles would help to guide the design and the opti¬ 
mization of artificial rheotactic active particles. 
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